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Abstract 

We discuss the magnetic phases of the Hubbard model for the honeycomb lattice both in two and 
three spatial dimensions. A ground state phase diagram is obtained depending on the interaction 
strength U and electronic density n. We find a first order phase transition between ferromagnetic 
regions where the spin is maximally polarized (Nagaoka ferromagnetism) and regions with smaller 
magnetization (weak ferromagnetism). When taking into account the possibility of spiral states, 
we find that the lowest critical U is obtained for an ordering momentum different from zero. 
The evolution of the ordering momentum with doping is discussed. The magnetic excitations 
(spin waves) in the antiferromagnetic insulating phase are calculated from the random-phase- 
approximation for the spin susceptibility. We also compute the spin fiuctuation correction to 
the mean field magnetization by virtual emission/absorpion of spin waves. In the large U limit, 
the renormalized magnetization agrees qualitatively with the Holstein-Primakoff theory of the 
Heisenberg antiferromagnet, although the latter approach produces a larger renormalization. 

PACS numbers: 71.10.Fd, 75.10.Lp, 75.30.Ds, 75.30.Kz, 81.05.Uw 
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I. INTRODUCTION 



The interest in strongly correlated systems in frustrated lattices has increased recently 
because of the possible realization of exotic magnetic states [1] , spin and charge separation 
in two dimensions [2], and the discovery of superconductivity in Naa;Co02.yH20 [3]. Many 
researchers have discussed superconductivity in non-Bravais lattices, mainly using self con- 
sistent spin fluctuation approaches to the problem [4-6]. The honeycomb lattice, which is 
made of two inter-penetrating triangular lattices, has received special attention after the 
discovery of superconductivity in MgB2 [7]. Additionally, the honeycomb lattice has been 
shown to stage many different types of exotic physical behaviors in magnetism and the grow- 
ing experimental evidence of non-Fermi liquid behavior in graphite has led to the study of 
electron-electron correlations and quasi-particle lifetimes in graphite [8]. 

Around a decade ago, Sorella and Tossatti [9] found that the Hubbard model in the 
half-filled honeycomb lattice would exhibit a Mott-Hubbard transition at finite U. Their 
Monte Carlo results were confirmed by variational approaches and reproduced by other 
authors [10, 11]. As important as the existence of the Mott-Hubbard transition in strongly 
correlated electron systems is the possible reahzation of Nagaoka ferromagnetism. The 
triangular, the honeycomb and the Kagome lattices were studied, but a strong tendency 
for a Nagaoka type ground state was found only in non-bipartite lattices (triangular and 
Kagome) [12]. On the other hand, the effect of long range interactions in half filled sheets of 
graphite was considered from a mean field point of view, using an extended Hubbard model. 
A large region of the phase diagram having a charge density wave ground state was found 
[13]. More recently, the existence of a new magnetic excitation in paramagnetic graphite 
has been claimed [14], but its existence was reanalyzed by two of the present authors [15]. 

In this work the magnetic phases of the Hubbard model in the honeycomb lattice are 
studied. In addition to the two-dimensional problem we also address the three-dimensional 
system composed of stacked layers. The critical lines associated with instabilities of the 
paramagnetic phase are obtained in the U, n plane (interaction versus particle density) . 
Spiral spin phases are also considered. A ground state phase diagram containing ferro and 
antiferromagnetic order is obtained. Interestingly, we find ferromagnetic regions with fully 
polarized spin in the vicinity of regions with smaller magnetization. The transitions from 
one to the other are discontinuous. 
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We also address the calculation of the magnetic excitations (spin waves) in the half- 
filled antiferromagnetic honeycomb layer within the random-phase-approximation (RPA). 
It is known that the Hartree-Fock-RPA theory of the half-filled Hubbard model is correct 
in both weak and strongly interacting limits: at strong coupling, the spin wave dispersion 
obtained in RPA agrees with the Holstein-Primakoff theory for the Heisenberg model; at 
intermediate interactions {U/t ~ 6), the RPA dispersion shows excellent agreement with 
experiment [16, 17]. The Hartree-Fock-RPA theory should, therefore, be considered as a 
usefuU starting point to study the intermediate coupling regime. Starting from the spin 
wave spectrum obtained in RPA theory, we calculate the quantum fiuctuations correction to 
the ground state magnetization arising from virtual emission/reabsorption of spin waves. In 
the strong coupling limit, we find a ground state magnetization which is about 67% of full 
polarization. This is not so great a reduction as predicted by the Holstein-Primakoff theory 
of the Heisenberg model, which is about 48%. 

Our paper is organized as follows: in section II we introduce the Hamiltonian and its 
mean field treatment. In section HI, we discuss the possibility of a well defined magnetic 
excitation in the paramagnetic phase. In the ordered phase at half filling, the spin wave 
spectrum is computed and the effect of different hopping terms in the spin wave spectrum 
is discussed. In section IV, the magnetic instability lines are obtained and the possibility of 
spiral spin phases for n < 1 is discussed. The corresponding lowest critical U is determined 
as function of the ordering wave- vector q. Section V is devoted to the phase diagram of 
the system, where two different types of ferromagnetism are found. The first order critical 
lines separating the three ordered phases are determined. Section VI contains a study of 
the renormalization of the electron's spectral function and magnetization by the spin wave 
excitations. 

II. MODEL HAMILTONIAN 

The magnetic properties of the honeycomb lattice is discussed in the context of the 
Hubbard model, which is defined as 

i,j,a i i,a 
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where Uj are hopping integrals, U is the onsite repulsion and // denotes the chemical poten- 
cial. The honeycomb lattice is not a Bravais lattice since there are two atoms per unit cell. 
Therefore, it is convenient to define two sublattices, A and B, as shown in Figure 1. 






^\ K 


r 


^ r 








Brillouin 



FIG. 1: Primitive vectors for the honeycomb lattice and the corresponding Brillouin zone. 

The expressions for the lattice vectors are 

ai = |(3,V3,0), a2 = |(3,-v^,0), a3 = c(0,0,l). (2) 

where a is the length of the hexagon side and c is the interlayer distance. The reciprocal 

lattice vectors are given by 

27r 27r In 

bi = — (1,^3,0), b2 = — (1,-^3,0), b3 = — (0,0,1). (3) 
6a 6a c 

The nearest neighbors of an atom belonging to the A sublattice are: 

5i = ^(l,V3,0) 52 = ^(1,-^/3,0) 53 = -ax 5" = ±cz (4) 

while the second nearest neighbors (in the plane) are: 5'^ — ±ai, 5'^ — ±a2, 5'^ = ±(a2 — ai). 
In a broken symmetry state, antiferromagnetic (AF) order is described by the average lattice 
site occupation: 

2 2 [-,jes 

where the axis ordering vector Q = (0, 0, Qz) will be used when studying multi-layers, n 
denotes the electron density, m is the staggered magnetization, and cr = ±1. We introduce 
field operators for each sublattice satisfying the usual Fourier transformations: 

fit - Jkihf.\ l\ - Jkiun (a\ 



(where N denotes the number of unit cells). Within a Hartree-Fock decouphng of the 
Hubbard interaction in (1) we obtain an effective Hamiltonian matrix 



Hii Hi2 




Oka 


H21 H22 




ha 



(7) 



ka 

with matrix elements given by 

i/n ^D{k) + U ^ , H^2 ^(t>k^ ^21, H22 = D{k) + U ^ 

where 



(8) 



<t>k^-tY, e'*'", D{k) = <t>'j, - 2t" cos{cK) - 1^, 0fe = -t' J2 e"'"' ■ (9) 

<5 <5' 

In the above equations t and t' are the first and second neighbor hopping integrals, re- 
spectively, while t" describes interlayer hopping. The dispersion relation for the case where 
t' = = is 



|0fe| = i-y^ 3 + 2 cos{VSaky) + 4 cos{3akx/2) cos{V3aky/2). 



(10) 



Diagonalization of the effective Hamiltonian yields a two band spectrum. The band energies 
are: 

TI I /TTm \ 2 

(11) 



U I /Um \ 2 



Because there are two sublattices, the Matsubara Green's function is a 2 x 2 matrix whose 
elements are given by: 



iuj ~ Ej (k) 



A . R* 



j^Ju-E^{k) 



where the coherence factors are: 



l^.,±(k) 



Uma ' 
2EJXj 



iuj — Ej (k) 



\Ba,±i^) 



(12) 
(13) 
(14) 
(15) 



1 + 



Uma 
2^±(k) 



2E±(k) 



(16) 
(17) 
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In the ferromagnetic (F) phase, the site occupation is the same for both sublattices: 

Tl Tfl 

< >= o + 77^ j^AB. (18) 



In this case the quasiparticle energy bands are given by 



U 



Elik) = D{k) + -(n - am) ± |0fc| 



(19) 



In the paramagnetic phase of the system the energies and propagators are simply obtained 
by setting m = in the equations above. The density of states of single electrons is shown 
in Figure 2 against particle density and energy. In the two upper panels we have included 
a second-neighbor hopping while in the two lower panels only nearest neighbor coupling is 
considered. An important feature is that p(e) vanishes linearly with e as we approach the 
half filled hmit, both for t' — Q and t' ^ 0. This is related to the X-points of the Brillouin 
Zone (see Figure 1), where the electron dispersion becomes linear: 

(dk denotes the deviation from the X-point). This dispersion is called the "Dirac cone". 




-3-2-101230 
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FIG. 2: Single particle density of states, p(e), for independent electrons in an honeycomb lattice. 
The left and right panels show p{e) as function of energy and electron density, respectively. The 
sohd line refers to t' = —0.2 and the dashed line to t' = 0. 



6 



III. COLLECTIVE EXCITATIONS AT HALF FILLING 



The magnetic excitations are obtained from the poles of the transverse spin susceptibihty 
tensor, x, which is definded, in Matsubara form, as 



xX'-{q, ^^n) = / dTe^^-'{TrSt{q, r)Sr{-q, 0)) 



(20) 



where i,j = a,b label the two sublattices (not lattice points) and S^{q), Sj'{q) denote the 
spin-raising and lowering operators for each sublattice. 

In the paramagnetic, F, or AF phases, the zero order susceptibility is just a simple bubble 
diagram with the Green's functions given in equations (12)- (15): 



ILOn — ilJJn 



(21) 



k,UJn 



Going beyond mean-field, the random-phase- approximation (RPA) result for the suscepti- 
bility tensor is obtained from the Dyson equation 

-1 



X = X° + Ux'^X ^ X 



X 



(22) 



where / denotes the 2x2 identity matrix. The poles of the susceptibility tensor, corre- 
sponding to the magnetic excitations, are then obtained from the condition: 



Det 



= 0. 



(23) 



We note that the tensorial nature of the spin susceptibility is a consequence of there being 
two sites per unit cell and is not related to the magnetic order in the system. 



A. Magnetic excitations in a single paramagnetic layer 

Here we discuss the possibility of existence of magnetic excitations in a single honeycomb 
paramagnetic layer. Our interest in this problem stems from a recent claim, by Baskaran and 
Jafari [14] , who recently proposed the existence of a neutral spin collective mode in graplicne 
sheets. In the calculations of Ref. [14] a half-filled Hubbard model in the honeycomb lattice 
(with t' — t" — 0) was considered but the tensorial character of the susceptibility was 
neglected [15]. Since inelastic neutron scattering can be used to study this spin collective 
mode in graphite, we decided to re-examine this problem taking into account the tensorial 
nature of the transverse spin susceptibility. 
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Collective magnetic modes with frequency cu and momentum q are determined from the 
condition (23) after performing the analytic continuation icu ^ cu + iO"^. The determinant 
is given by 



D^_{q, a;) = 1 - 2Uxfr + [{x^-r ' xV-''xT\ , (24) 

where we have taken into account that in a paramagnetic system x^l"" = X+-^- Below the 
particle-hole continuum of excitations, the spectral (delta-function contributions) part in 
X+-''(^;'^ _l_ vanishes and there is the additional relation = {x^+-^)* ■ Collective 

modes are only well defined outside the particle-hole continuum (inside the continuum they 
become Landau damped). Wc scarchcd[15] for well defined magnetic modes, oj{q)^ below the 
continuum of particle-hole excitations, and found no solutions for any value of the interaction 
U. In Figure 1 of Ref. [15] we plot D^_{q,uj) for eight different q- vectors and u ranging 
from zero to the point where the particle-hole continuum begins. Our analysis reveals that 
the full tensorial structure of the Hubbard model's RPA susceptibility in the honeycomb 
lattice docs not predict a collective magnetic mode. 



{0)aa\2 (0)a6 (0)6a 



B. Spin waves in the antiferromagnetic layer 

The spin wave dispersion uj{q) for the AF layer with one electron per site can be obtained 
from equations (21) and (23) using expressions (12)-(15) for the propagators. Spin wave 
spectra, for different values of second- neighbor hopping, t' , are ploted in Figures 3 and 4. In 
the large U limit, spin wave energies agree with those obtained from the Holstein-Primakoff 
theory of the Heisenberg model. We give an analytical derivation of this limit in Appendix 
B. The Holstein-Primakoff result for the Heisenberg model in the honeycomb lattice, which 
is derived in Appendix C, can be written as 

uJHp{q) = JS^z^-\(t>{q)\'' . (25) 

This result can be mapped on the Hubbard model provided that J = At^/U and S = 1/2. 
Figure 3 shows the spin wave energies for the 2D lattice {t' = 0) along a closed path in the 
Brillouin Zone. Energies in Figure 3 are normalized by the Holstein- Primakov result at the 
K-point, ujhp{K) (see Figure 1). It can be seen that the results for [/ = 8 are very close 
to the asymptotic behavior of the RPA, whereas, for smaller [/, the spin wave energy is 
reduced. The effect of t' on uj{q) is depicted in Figure 4. It is of particular interest the fact 
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that the dispersion along the X — K direction is almost absent for [/ > 4. The presence of 
t' does not change this effect. 




r X K 



FIG. 3: Spin-wave excitation spectrum for several values of U . The dashed-dotted line gives the 
Holstein-Primakoff result for the Heisenberg antiferromagnet in the honeycomb lattice. 
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FIG. 4: Spin-wave excitation spectrum for several values of U and ^ 0. 



IV. MAGNETIC INSTABILITIES 



The magnetic instabilities in the paramagnetic phase can be obtained from the divergence 
of the RPA susceptibihties at critical values of the interaction, [4, driving the system towards 

9 



a magnetically ordered phase. At a given electron density n we always find two instability 
solutions, one ferromagnetic and one antiferromagnetic. One of these solutions minimizes 
the free energy. Since Uc is determined from D+_(q,0) = 0, taking into account that 

and x^l"^ = {x+-'^)* the paramagnetic phase, we obtain: 

1 



(0)aa 
X+- 



K+- 



^'^ (0)oa I I (0)ab , 



(26) 



Figure 5 shows Uc obtained from the static uniform susceptibilities {q — and u; — 0), as 
a function of electron density for various values of f. Detailed equations for the instability 
lines are given in Appendix A. The left panel of Figure 5 refers to the 2D case, corresponding 
to a single honeycomb layer, whereas the right panel refers to the 3D system with a constant 
interlayer hopping t" = 0.1. The Van- Hove singularity (associated with the X point) plays 
an important ] 




1.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0. 
n 



1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
n 



FIG. 5: left panel: Effect of t' on the instability lines, as determined from the equation (26), 
for a single honeycomb layer, right panel: Effect of t' on the instability lines, as determined 

from equation (26), for a layered honeycomb. This panel differs from the other inasmuch a small 
t" =0.1 hoping term was included coupling the 2D layers. 

As we have already mentioned, the two solutions of Eq. (26) correspond to two different 
magnetic transitions, one between a paramagnetic phase and a ferromagnetic phase and 
another between a paramagnetic phase and an antiferromagnetic phase. That this is so can 
easily be confirmed by solving the self-consistent equations for the ferromagnetic and the 
antiferromagnetic magnetizations, respectively, derived from the HF Hamiltonian (7). By 
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minimizing the free energy with respect to magnetization, one finds the following expressions 
for ferro and antiferromagnetic magnetizations 

where f{x) is the Fermi function and Cfe — U'mAF/{'^\(pk\)- Letting both mp and ttiaf 
approach zero, one obtains the same lines as those in Figure 5. Generally speaking, for 
electron densities lower than 0.85, the value of Uc that saparates the paramagnetic region 
from the ferromagnetic region is lower than the corresponding value of Uc separating the 
paramagnetic region from the antiferromagnetic region. The critical U associated with the 
ferromagnetic instability increases with t'. The size of the paramagnetic region in Figure 
5 increases with t'. On the other hand, for t' — 0.2, we see that the critical fine for the 
ferromagnetic region is very close the critical line of the antiferromagnetic region. Therefore, 
the ferromagnetic region is progressively shrinking with increasing t'. If we now turn to 
densities larger than 0.85, we find that the antiferromagnetic critical line is the one with 
lowest Uc- However, in contrast to lower densities, the antiferromagnetic critical line hardly 
changes when varying t'. This description applies equally well to the single honeycomb layer 
and weakly coupled layers, even though the quantitative functional dependence of Uc on n 
is different in the two cases, the main difference coming from the van Hove singulary present 
in the 2D case. At finite temperature the van Hove singularity is rounded off and the 2D 
phase diagram will be much more similar to the 3D case. We therefore, consider that a weak 
3D inter-layer coupling does not qualitatively modify the conclusions valid for the 2D case. 

Besides coUinear spin phases, the system may also present non-coUinear - spiral - spin 
phases in some regions of the phase diagram. We now study what are the changes in the 
critical U values determining the instability of the paramagnetic phase if we allow for non- 
coUinear ground states, since it is well known that the Hubbard model on bipartite and 
non-bipartite lattices can have the lowest Uc for spiral spin phases [12, 18, 19] for some 
electronic densities. In a spiral state, the spin expectation value at site i, belonging to 
sublattice u = a,b, is given by [20] 

{S'^) = ^{cos{q-R^),sm{q.R-)). (28) 

If qr 7^ 0, the ferromagnetic and antiferromagnetic spin configurations become twisted. We 
shall refer to the twisted q configurations as 'F^' whenever = m^, and 'AF^' 
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whenever tua = —mB- The criterion for choosing the q- vectors is taken directly from the 
geometry of the lattice by requesting a constant angle between spins on neighboring sites, 
i.e. q • 6^ = q ■ §2 = —q • 6^. Unfortunately, however, this cannot be achieved in the 




FIG. 6: (color on line) Fq (upper) and AFq (lower) spin configurations for = f • 

honeycomb lattice with only one q-vector. The closest one can get to a 'true' spiraling 
state is by letting q ■ <5i = —q • 8^ (or equivalently, q ■ 62 = —q ■ S^), which implies that 
q = {qx, qy) = qx{^, :^) (q = qx{^, ^))- the moment we let qz be zero which means that 
we consider identical layers. The condition q ■ 5i = —q ■ 63 means that the increase in spin 
angle between two lattice sites in the — ^3 direction is the same as the increase in spin angle 
between two lattice sites in the direction. There is no increase in the spin angle in the 
62 direction. Examples of the spin-configurations obtained in this way are shown in Figures 
6 and 7. Several notes are in order at this stage. First, although we do not have a 'true' 
spiraling state over the whole lattice, we do have a spiraling configuration in the — ^3 and Si 
directions, as can be seen from the Figures 6 and 7, going from the lower left to upper right. 
Secondly, when travelling along the ^2 direction, the spin angles do not increase. Instead, 
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FIG. 7: (color online) Fq (upper) and AFq (lower) spin configurations for qx = 3-- 

neighboring spins in this direction are always aligned ferromagnetically when = ms-, and 
antiferromagnetically when itia = —ttib- However, two successive 62 bonds ('sliding down' 
the lattice from left to right) have the same increase in spin angle as any two neighbors 
connected by —S3 or 61. The g-vector (i.e. the spin configuration) that a system with a 
given density would prefer is the one with the lowest value of Uc{q)- In Figure 8 we present 
a curve showing the q vectors that minimize Uc{q), as functions of particle density n. We 
consider discrete values = with i = 0,1, 12. The dependence on t' is overall the 
same as that discussed for q = (for example, the shrinking effect with increasing t' is 
also seen here). There is no reason to restrict q to integer multiples of j^, other than a 
pure computational one. By performing the same calculation with more g-vectors, the 'step 
function' like appearance of the lower graphs of Figure 8 can be smoothed out. Our analysis 
is sufficient, however, to get an insight into how the q vectors (which minimize Uc) vary with 
n. 

The solid line limiting the paramagnetic region is shown in the lower graphs of Figure 
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Minimum RPA Lines 2D Honeycomb 



Minimum RPA Lines 3D Honeycomb 




FIG. 8: The upper panels show the minimum Uc{q) according to Eq. (26). The sohd hne separates 
the paramagnetic phase from magnetically ordered phases, while the dashed line separates differ- 
ently ordered magnetic phases. The lower panels show the Qx component of the ordering vector 
(corresponding to the minimum Uc) as function of electron density n. 

8). We see that the behavior of q^, as function of n, is almost the same for the 2D and 3D 
cases. As the system approaches half filling, the prefered spin configuration approaches that 
with q = 0. In a doped system, however, minimization of Uc{q) is attained for a non-zero 
q. It is also seen that the dependence of q on n is not monotonic. Either in 2D or 3D, 
goes all the way from (at n = 1) to tt, displaying two local maxima (and a local minimum 
in between) as n ranges from 1 towards 0. 

The value of q^ reaches a local minimum a,t q^ = at n — 0.37 (in 2D) or at n = 0.45 
(in 3D). For even lower densities, q^ attains another maximum at q^ = vr, which means 
that the spins of any two nearest neighbors, in the —ds and <5i directions, point exactly in 
opposite directions to each other. 

The same type of behavior is seen also for the critical line separating magnetically ordered 
phases (dashed hne). Again, the 2D and the 3D cases are very similar to each other. For 
densities around 0.30 — 0.35 (2D) and 0.35 — 0.40 (3D), we have qx — ^ yielding the lowest 
Uc- Moreover, the solid and the dashed lines coincide, illustrating the previously mentioned 
ferromagnetic 'shrinking out' effect. In other words, for qx — the two solutions of 
Uc{q) almost coincide for all n, leaving only a thin strip of ferromagnetism between the 
paramagnetic and the antiferromagnetic regions. Although this is true for all n, it is only 
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for n = 0.37 - 0.40 (2D case) and n = 0.35 - 0.37 (3D case) that Uc{qx = f|) is minimum. 

So far, our analysis has been restricted to qr-vectors lying in the x — y spin plane. This 
means that two inter-layer neighbors have the same spin. If we now consider neighboring 
layers with opposite spin, we put = tt. At half filling, the lowest Uc{0, 0, tt) = 2.04 limiting 
the paramagnetic region is lower than the corresponding Uc{0, 0, 0) = 2.35, independently of 
t'. Moreover, for n = 1, Uc{qx, ^j'^) always lower than Uc{qx, ;^)0) for Qx, showing 
that, at half filling, we should expect antiferromagnetic ordering along the 2;-direction. 

The study above was focused on the second order instabihty fines, both in the case of 
coUinear and spiral spin phases, being clear that spiral states have a lower critical U— value, 
over a large range electronic densities. It is instructive to compare our results with those 
of Ref. [12]. Looking at Fig. 2 of Ref.[12] we see that for the triangular lattice there are 
some finite regions where the more stable ground states correspond to spiral states. These 
regions are located at electronic densities smaller than 0.5 and larger than 0.8. Since the 
honeycomb lattice consists of two inter-penetrating triangular lattices we expect the same 
type behavior, at least at the qualitative level. That is, we do expect to have finite regions of 
the phase diagram where spiral phases have the lowest energy. Also, in Ref. [12] the authors 
do not discuss the full phase diagram of the Hubbard model in the honeycomb lattice, as we 
do in next section. They are primarily interested in the stability of the Nagaoka state. Their 
study is done using three different approaches (i) The Hartree single fiip ansatz; (ii) the SKA 
Gutwiller ansatz; (iii) the Basile-Elser ansatz. A comparison can be established between the 
the Hartree single fiip ansatz which roughly speaking, produces a straight line for all densities 
at the on-site Coulomb interaction [/ ~ 5, and our self consistent Hartree- Fock study. If we 
forget, for a moment, the van Hove singularity, both results are qualitatively the same for n 
up to 0.8. Above this value our Hartree-Fock analysis, forgetting about the existence of the 
antiferromagnetic phase, predicts a very strong increase of the critical U value (not shown in 
Fig. 5, since the AF phase presents the lowest critical [/-value), in agreement with the SKA 
ansatz. This behavior is not captured by the the Hartree single fiip ansatz. It seems that 
our study interpolates between the Hartree single fiip ansatz for low densities and the SKA 
ansatz for densities above 0.8. Quantitatively there are differences between the two studies, 
which are understandable on the basis of the different types of proposed ground states. 
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V. PHASE DIAGRAM 



As we mentioned in the previous section, the study of Ref. [12] is mainly concerned with 
the stabihty of the Nagaoka state, and in the previous section we studied the values of the 
Hubbard interaction associated with instabilities of the paramagnetic system. The transition 
from the paramagnetic to a magnetically ordered state is determined by the lowest Uc- Since 
we have found the possibility of having, at least, two (ferro and antiferro) different types of 
ground states, then in the case where interaction is stronger than both critical values, we 
need to address the problem of competition between the two ordered phases. The phase with 
the lowest free energy is the one prefered by the system. In this section we restrict ourselves 
to the study of a single layer but we shall consider different band structures. Spiral states 
will not be considered, since we are most interested in a weak ferromagnetic phase showing 
up in region of the phase diagram where the studies of Ref. [12] suggest that the collinear 
ferromagnetic (fully polarized) phase should be the most stable one. In the ferromagnetic 
phase we distinguished two types of ferromagnetic ground states: the Nagaoka ground state, 
with a maximally polarized spin {mp — n), and a weak ferromagnetic state with mp < n. 
The order parameter and free energies were obtained from the mean field Hamiltonian (7). 
Figure 9 shows the § 
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FIG. 9: left panel: Ground state phase diagram of the Hubbard model in the (n, U) plane for a 
single layer with t' = 0. right panel: Ground state phase diagram of the Hubbard model in the 
(n, U) plane for a single layer with t' = —0.2. In both cases dashed and continuous lines represent 
first and second order transitions, respectively. 



16 



The effect of t' on the phase diagram can be seen in right panel of Fig. 9. In Figure 9 
the dashed hnes represent first-order phase transitions, where the order parameter do not 
vanish smoothly, while continuous lines represent second order transitions, where the order 
parameter vanishes smoothly, but its first derivative is discontinuous. In both cases {t' = 
and t' 7^ 0) we find a finite region of weak ferromagnetism. In general the Nagaoka phase is 
more stable for large U. The weak ferromagnetic phase is separated from the Nagaoka phase 
by first or second order transition lines, depending on the path followed on ([/, n) diagram. 
The second order transition manifests itself through a discontinuity of the derivative of the 
magnetization with respect to U. At n = 0.75 the instability line towards the ferromagnetic 
phase shows a dip (pronouced if t' = 0), which is due to the logarithmic van- Hove singularity 
at n = 0.75. A negative t' produces two effects on the phase diagram: (i) the instability line 
towards the F phase moves downwards; (ii) the point where the instability lines towards 
F and AF meet moves to larger n. Similarly to what was found in the previous section, 
the overall effect of t' is to modify the ferromagnetic region of the phase diagram. Further, 
for negative t' we expect coUinear ferromagnetism to exist over a large phase of the phase 
diagram relatively to the case t' > 0, since it is well known that a negative t' stabilizes the 
ferromagnetic phase. On the other hand we don't expect the phase diagram presented in 
this section to be fully accurate for low densities, where the findings of Ref. 12 should apply. 

The first order critical lines do separate two different ferromagnetic (or ferromagnetic from 
antiferromagnetic) regions, in what concerns the total magnetization. In view of the results 
published in Ref. [21], where a first order transition between the two competing phases 
is transformed by disorder into two second order phase transitions, we expect the same 
behavior to apply here, that is, disorder may change the order of the transition, since the 
arguments put forward in Ref. [21] are of very general nature. It would be very interesting 
to study whether the introduction of disorder in the system could change the nature of the 
first order transitions. 

VI. QUANTUM FLUCTUATIONS 

This section is devoted to the calculation of quantum fluctuation corrections to the mag- 
netization. An analogous calculation for the Hubbard model in the square lattice in the 
t/U ^ limit was skecthed by Singh and Tesanovic.[22] 
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The computation of the renormahzed staggered magnetization requires the evaluation 
of the Feynman diagram shown in Figure (10), which shows the second order (in the in- 
teraction U) contribution to the self-energy. The diagram describes the emission and later 

4 




Pf p-q| Pf 



FIG. 10: The self-energy for a |-spin electron. The bubble represents the transverse susceptibility 
computed in RPA. 

absorption of a spin wave by an up-spin electron. The emission and absorption processes are 
accompanied by electron spin reversal. This effect, consisting of virtual spin flips, is going 
to renormalize the staggered magnetization. The spin-| electron Green's function is 

G^{p, iuj) = ^°(p, iuj) + ^°(p, it^)S|(p, iuj)Q^{p, iu) , 

hence, G^^ = [G^'^^]^^ — S^^. Here, G^ denotes the Hartree-Fock Green's functions matrix 
appearing in equations (12)-(15). The self-energy matrix is given by 

S|^ (p, zu;) = U'^ Yl (P - - ^f^)x^T^^'^ (q, , (29) 

where i,j are sublattice indices. The self-energy for a J,-spin electron would be similar to that 
in equation (29) with the ^^^-'^-'-spin reversed and X-+ repaced with x+~- The renormahzed 
staggered magnetization at T = is given by 

^ = -^E r ^^[I^G::Uk,u:)-ImGl%Ak,uj)], (30) 

where ImG^^-^^^^{k,uj) stands for the imaginay part of the retarded Green's function for a 
spin (J electron. 

The RPA susceptibility has poles corresponding to the spin waves calculated in section 
III, with energy ^ |0(k)p/[/, but it also has poles describing a particle-hole continuum 
of excitations at higher energies (of order U). In what follows we ignore this particle- hole 
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continuum and take into account only the contribution from the spin wave poles to the 
selfenergy. Physically, this means that we shall calculate the magnetization renormalized by 
the spin waves. To this end, we start by replacing the susceptibility in equation (29) by the 
expression 

= + (31) 

luj — a;(q) luj + a;(q) 

where R'^[±.uo{c\ j\ denotes the residue of x^-+^^^^ at the spin wave pole with dispersion a;(q). 
Equation (31) describes an effective spin wave propagator. After performing the Matsubara 
frequency summation in equation (29) we obtain: 

. ,_U^x- \ "^"^{g{°-^'(P - q)}^^^[-^(q)] num{g{°y(p - q)}it:-^[a;(q)) ^^^ 
T ^P' ''^) - N ^ I iu; + uiq) + E+(p - q) iu - u;(q) - E+(p - q) ^l ' 

where we have introduced the notation num{^^°^'-' } for the numerators of the Green's func- 
tions, as expressed in equations (12) through (17). 




FIG. 11: The magnetization in the half-filled honeycomb AF layer. The continuous line represents 
the Hartree-Fock result. Renormalized magnetizations are shown for different lattice sizes: 20x20; 
62x62; 82x82. The vertical dashed line represents the mean field critical U value at which the 
magnetic instability develops. 

Figure 11 we show the renormalized magnetization versus U. The Hartree-Fock magne- 
tization is also shown in the Figure 11 for comparison. The calculation was performed for 
three different lattice sizes. It can be seen that convergence does not require a very large 
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number of k points in the Brillouin Zone. This is not surprising because the Hartree-Fock 
magnetization itself aheady converges to the correct value in a 40x40 lattice. We have 
also checked that the RPA propagators return the original electron density n = 1, meaning 
that no spectral weight was lost in the used approximation for the self energy. In the large 
U limit, the renormalized magnetization saturates at about 67% of the (fully polarized) 
mean field value. This is in qualitative agreement with the Holstein-Primakoff result for the 
S — 1/2 Heisenberg model in the honeycomb lattice, which predicts a ground state magne- 
tization of 48%. We should remark, however, that the spin wave spectrum calculated within 
RPA theory has shown much better agreement with experimental results for Mott-Hubbard 
antiferromagnctic insulators than the Holstein-Primakoff theory [16, 17]. 

In Figure 12 we show the imaginary part of the electron's Green's function at negative 
frequencies, on both sublattices, for two different values of U. It is clear that, for strong 
couphngs, part of the Hatree-Fock spectral weight is shifted to the bottom of the (negative) 
energy band. This shifting of the spectral weight is responsible for the renormalization of 
the staggered magnetization. It is interesting to see that for low U the spectral weight is 
most significant at high energy, in the interval [-2,0 [, with a much smaller weight in the 
interval ]-4,-2[. At a stronger Hubbard interaction most of the high energy spectral weight 
(previously in the interval [-2,0[) has been displaced to lower energies and become localized 
around well defined energies, whereas the spectral weight at intermediate energy (in the 
interval ]-4,-2[) remains essentially unchanged. Therefore, increasing Hubbard coupling has 
the efect of displacing the distribution of spectral weight from the top to the bottom of the 
energy band. 

Finally, a comment regarding approximation (31). The commutation relation between 
the spin raising and lowering operators, 

p,p' p 

is equivalent to the following relation between the Hartree-Fock magnetization, m, and the 
transverse susceptibilities: 

X»A(q, r = 0+) - x"A(q, r ^ 0') ^ j> -^X^-We'^' - f -^X%{z)e^'" 

— —m , (33) 

at T = 0. The integration of the term e"^^'^ (e^°^) is performed along the semi-circular 
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FIG. 12: Imaginary part of the retarded electron Green's function multipied by —1, —ImG^"'^^^^ (w), 
versus negative frequency. The Green's function includes the quantum fluctuations. 



contour on the right (left) half of complex plane. Approximation (31) would predict 



(34) 



Indeed, we have checked that our numerical calculation of the residues satisfies (34) to an 
accuracy of 1.3%. 

VII. FINAL REMARKS 

In this paper we have studied the magnetic properties of the Hubbard model in honeycomb 
layers. Our study focused on the instabilities of the paramagnetic phase, on the magnetic 
phase diagram and on the collective excitations of the half filled phase. Of particular interest 
is the fact that it is not possible to describe a true spiraling state in the honeycomb lattice, as 
opposed to the usual cubic case. As a consequence, the magnetic spiral order follows a kind 
of one dimensional path over the 2D lattice. This kind of ordering, here studied at mean 
field level, may have important consequences to the study of spin charge separation in 2D 



21 



lattices. Also interesting, was the identification of two types of ferromagnetic order, which 
have eluded previous studies. For moderate values of U and electron densities not far from 
the half filled case, a region of weak ferromagnetism was found to have lower energy than 
the more usual Nagaoka ferromagnetic phase. The renormalization effect of the spin wave 
excitations on the Hartree Fock magnetization was also studied. However, our calculation 
does not take into account the renormalization of the mean field critical U. It is well known 
that quantum fiuctuations should induce an increase the value of Uc- Our calculation cannot 
capture this effect, since it only takes into account the effect of well defined spin waves. We 
believe, however, that the calculation can be extended to include the effect of high-energy 
damped particle-hole processes leading to a renormalization of Uc, but this would require a 
modification of our numerical calculations and a significant increase of the computational 
time. 



APPENDIX A: USEFUL EXPRESSIONS FOR THE Uc CRITICAL LINES AT 



In this appendix, we derive the equations for the critical lines from the static suscepti- 
bilities (q = and u = 0). Our starting point is the zero order spin-spin susceptibility in 
equation (21). The Green's functions in the paramagnetic region are obtained from equations 
(12)-(15) after setting the magnetization to zero. Performing the Matsubara summations in 
(21), the analytical continuation and taking the zero frequency limit, we obtain 



xK(Q,0) = ^^(M++(fc,q) + M+_(fc,q) + M_+(fc,g) + M__(fe,g)) (Al) 

k 

xK(9,0) = ^^e^(^'=--^'=)(M++(fc,qr) -M+_(fc,g) -M_+(fc,q) + M__(fc,q)](A2) 
fc 

where i/jk — arg(0fe). The critical interaction strength, Uc , is given by Uc/N — [x+-fi ± 
|x+^"ol ™ limit qr — > 0. Expanding all q dependent quantities around the point 
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q = up to first order, we obtain 
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Inserting this result in the expression for Uc gives (for q — 0): 



(§) =lj2{^[E4k)]+S[E4k)]} (A6) 



l^ g(|0fe|-|D(fc)|) 



'?7eV' _ ^yO{\(l>k\-\Dik)\) 

k 

We recognize the density of states, p(e) = jj: ^k{^{E+{k) + // — e) + S{E_{k) + fji — e)}, 
appearing in equation (A6), which is just the Stoner criterion. The critical interaction 
strengths are given by 

2 

iV Z^fe |<;!,fe| 

Note that all t' and dependence is contained in D{k). Of course, these equations could 
also have been obtained by taking the limit mp, tuaf — > in equation (27). 



APPENDIX B: LARGE U RESULTS FOR THE SUSCEPTIBILITIES AND SPIN 
WAVES 

We give asymptotic expressions for the susceptibilities x^j._ (z, q) and spin wave dis- 
persion for a half-filled honeycomb antiferromangetic layer with nearest neighbor hop- 
ping. In this case, the chemical potential = and the two energy bands are given by 

^(k)± = ±^(^)V|0fcP. 

The expressions for coherence factors appearing in the single electron propagators, ex- 
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panded up to second order in t/U , are: 



iA.+(k)r 




(k)r 


= isT,-(k)|2« 


IA,-(k)P 


= IA,,+(k)p = IB,,_ 


(k)r 


= |5T,+(k)r«^ 


^I,_(k)5,,_(k) 


= -^;,^(k)5T,+(k) 












,(k)Si,+(k) 



1 - 



^*(k) 

We therefore may use the aproximate expressions for the x° _ susceptibihties: 

|0(k)p + |0(k + q)pN 



(Bl) 
(B2) 

(B3) 





--E 

k 




^ k ^ 


X(°)''»(^,q) «^ 






k 



1 - 



(k)p + |0(k + q)p 

1 



(B4) 



) 



(B5) 



m 0*(k + q 



Z-E(k) -E(k + q) ;2 + E(k) +E(k + q)/ f/W 

1 X0*(k)0(k + q 



B6) 



-B7) 



.;2-E(k) -E(k + q) ^ + E(k) +E(k + q)/ ^72^2 

We anticipate that the spin wave energies are of order z /U so that we may use the 
expansion 



1 



z + ^(k) + ^(k + q) Um 



1 



(k + q) 



Um 



+ ... 



in equations (B4)-(B7). The condition (23) for the spin wave dispersion now takes the form: 



C/2m4 



1 + 



C/2m6 I 



-5]0*(p)0(p + q) 



(B8) 



p p 

But we must take into account that the self-consistent equation for the Hatree-Fock magne- 
tization, expanded to second order in t/U, is 

2 



i-i 



(B9) 



Introducing (B9) in (B8) we finally obtain the spin wave dispersion: 

' ^ ^^^^ " E i<^(p) if - 1 ^ E <^*(p)'^(p + r ' (BIO) 

V p p 

which agrees with the result predicted by the Holstein-Primakoff theory. 
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APPENDIX C: HOLSTEIN-PRIMAKOFF ANALYSIS OF THE HEISENBERG 
MODEL 

The Heisenberg Hamiltonian in the honeycomb lattice is given by 
^ = ^ E ^^^^Us + liS^ST^, + S.St,,)] + ^ E ^^i^Us + liS^ST^s + S^Sts)] ■ (CI) 

ieA,5 ieB,S 

We introduce two sets of operators 



S^ = -alui + S, S+ = \j2S-a\ai , S+ = a\ \l2S-a\ai , (C2) 

and 

St^-blk + S, St^^2S-b%k, St^h\^j2S-h\h. (C3) 

Making the usual linear expansion and introducing the momentum representation for the 
bosonic operators, the Hamiltonian can be written as 

H = -JNazS^ + JzS E(«fc«fc + blbk) + JS J2i^ik)akb-k + {k)bl^ai) . (C4) 

k k 

Next we introduce a set of quasiparticle operators defined by 

«fe = ^kllk - ^kl2,k , blf^ = Uk^lf^ - vl'yi^k , (C5) 

where the coherence factors obey \uk\'^ — \vk\'^ ~ 1. After introducing the above transforma- 
tions in the Hamiltonian we find 

H = -JNazS^ + J2('^Jz^\'^k\'^ - J^^l^i^WK- JS(l)*(k)vluk) 

k 

+ E [J^S{\uk\^ + \vk\^) - JS(P{k)vkul - JS(P*{k)vluk)hlkli,k 

k;i=l,2 

+ Y^[{-2JzSvkUk + JS(t){k)vkVk + JS(t)*{k)ukUkWi^kll,k + H.c] , (C6) 

k 

which implies the conditions 

JzS{\uk\'^ + \vk\'^) - JS(f){k)vkul - JS(j)*{k)vluk) = uj{k) , 

—2JzSvkUk, + JS4>{k)vkVk + JS(p*{k)ukUk = 0. (C7) 

The second condition reveals that we can choose Uk to be real and Vj^ — 0(fc)Q;(fe), with 
a{k) real. After some straightforward manipulations we find 

1 ~ TO 

.(k) = JSV^^^, = -_ + __. (C8) 
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The stagered magnetization is given by 

znB[uj{k)] 



(C9) 

and at zero temperature we assume nB\^{k^\ = 0. Computing the integral gives a magneti- 
zation value of 0.24, that is about 50% the Neel value |. 
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